Positional Order and Diffusion Processes in Particle Systems 

Hiroshi Watanabe,^* Satoshi Yukawa,^ and Nobuyasu Ito^ 

^Department of Complex System,s Science, Graduate School of Information Science, 

Nagoya University, Furouchou, Chikusa-ku, Nagoya 464-8601, Japan 

^ Department of Earth and Space Science, Graduate School of Science, 

Osaka University, 1-5 Yamadaoka, Suita, Osaka 565-0871, Japan and 

•^ Department of Applied Physics, School of Engineering, 
The University of Tokyo, Hongo, Bunkyo-ku, Tokyo 113-8656, Japan 

Nonequilibrium behaviors of positional order are discussed based on diffusion processes in particle 
systems. With the cumulant expansion method up to the second order, we obtain a relation between 
the positional order parameter ^ and the mean square displacement M to be ^ ~ exp(— K^M/2d) 
with a reciprocal vector K and the dimension of the system d. On the basis of the relation, the 
behavior of positional order is predicted to be 'I' ~ exp(— K Dt) when the system involves normal 
diffusion with a diffusion constant D. We also find that a diffusion process with swapping positions 
of particles contributes to higher orders of the cumulants. The swapping diffusion allows particle to 
diffuse without destroying the positional order while the normal diffusion destroys it. 
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The melting behavior of the hard-disk system was re- 
ported first by Alder et al. Q, and they showed that 
only repulsive interactions can involve the melting tran- 
sition. This melting transition also confirmed in three- 
dimensional systems and is now often referred to Alder 
transition. However, Mermin ruled out the positional 
long range order in two-dimensional particle systems ^ . 
Therefore, the melting processes of two-dimensional sys- 
tems are different from that of three-dimensional sys- 
tems. Halperin, Nelson, and Young proposed the two- 
dimensional melting theory Q based on Kosterlitz Thou- 
less transition |4|, and Chui proposed another theory 
predicting the first order transition based on the grain 
boundaries excitation Q. While many researchers have 
been studying this problem ^ H la H, llfl Hll i the na- 
ture of two-dimensional melting has been still a matter 
of debate [Tjl- So far, most of numerical works have fo- 
cused on the equilibrium state of the system mainly using 
Monte Carlo methods. Recently, the nonequilibrium be- 
haviors of the bond-orientational order parameters has 
been studied to obtain the equilibrium properties of the 
hard-disk systems 13j. These studies are based on a new 
strategy for the simulation, called Non-equilibrium relax- 
ation (NER) method J^. Zahn and Maret studied time 
dependent parameters in two-dimensional colloidal parti- 
cle systems Tsj. They pointed out that static properties 
are not appropriate measure to distinguish between the 
solid and the fluid, since the mean square displacement 
diverges very slowly. Therefore, it is necessarily to study 
the dynamic behaviors of order parameters in the parti- 
cle systems. In this letter, we study the dynamics of the 
positional order parameter based on diffusion processes. 
We also treat two- and three-dimensional systems at the 
same time, since many studies have focused only on the 



* E-mail: hwatanabe@is.nagoya-u.ac.jp 



two-dimensional melting and, to our knowledge, there are 
less studies about the three-dimensional positional order. 
Consider a d-dimensional system with N particles. A 
positional order parameter ^ of the system is defined to 
be 



1 ^ 
* = -^expHKT,), 



(1) 



with the position of the particles r^ and one of the recip- 
rocal vectors K of the system. Let Rj be the equilibrium 
position of the particle i and u.^ the deviations from it, 
namely, r^ = R^ -I- u^. The positional parameter is re- 
duced to be, 



* = (exp(-iK • xi-j)) . 



(2) 



where (• • •) means the average for all particles. Assuming 
that the all components of u^ have the Gaussian distri- 
bution, Eq. |(2Jl is reduced to be Il7l|, 

vI/ = exp(-l/2((K.u,)2)). (3) 

Assuming that u^ is isotropic, we have, 

{{K.u.y) = K^u^)/d iK^\K\). (4) 

From Eqs. © and Q), we obtain the relation between 
the positional order and the diffusion to be, positional 
order parameter to be. 



^ 



exp 



K^u^y 



2d 



or equivalently, 



(u?) 



2d 



In*. 



(5) 



(6) 



Note that, the above argument is the cumulant expan- 
sion. The positional order parameter is the characteristic 
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FIG. 1: Time evolution of the mean square displacement 
M = (ui ) in (a) two- and (b) three-dimensional systems. The 
decimal logarithm are taken for the both axes. The dashed 
lines denote the diffusion in the low density limit p^ . Number 
of particles N = 23288 for two- and N = 36000 for three- 
dimensional system with the periodic boundary condition. 



FIG. 2: Time evolution of the positional order parameter 
and calculated values from the diffusion for (a) two- and (b) 
three-dimensional systems. The solid lines are the positional 
order parameters and the symbols are the calculated values 
using Eq. l|^. It shows good agreement in the region where 
the positional order parameters are not so small. 



function of displacements. Assuming the distribution of 
the displacement to be the Gaussian distribution, we can 
express the positional order parameter only with the sec- 
ond order cumulant, which is diffusion. 

When a system involves the normal diffusion, the 
asynptotic behavior of the mean square displacement is 
expected to be. 



(u^) - 2dDt, 

with a diffusion constant D. From Eqs. ll^Jl and iQ, 
asymptotic behavior of the positional order to be. 



(7) 
the 



(8) 



regardless of the dimension. It implies that when the 
system involves the normal diffusion, the positional or- 
der should decay exponentially with the decay time D~^. 
This limits the diffusion behavior in solid phases. In the 
solid phase of the system with d > 3, the parameter ^ 
has non-zero value in the equilibrium state. Therefore, 
the mean displacement cannot become larger than some 
constant value. The behavior in two-dimensional solid 
is different from those in d > 3. On the basis of the 



Halperin-Nelson- Young theory |3], the positional order 
parameter in two-dimensional solid behaves as, 



^{t) ^ t- 



(9) 



The mean square displacement in two-dimensional solid 
behaves logarithmically as, 



(u?) 



4A 
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(10) 



Therefore, two-dimensional solid cannot involve the nor- 
mal diffusion process in the usual sense. 

The above arguments are based on the cumulant ex- 
pansion up to the second order. In order to check the 
validity of our arguments, we perform numerical simula- 
tions. For the simplicity, we treat the hard-particle sys- 
tems. There are two kinds of methods to study time evo- 
lution of a particle system, a molecular dynamics (MD) 
method and a Monte Carlo (MC) method. The MD simu- 
lation is performed by integrating the classical equations 
of motion. In the hard particle system, the time evolution 
is performed by proceeding collision events. This algo- 
rithm is called the Event-Driven method, which is very 



efficient to treat the hard-particle system [T^ [T^ |20| . 
When the time evolution of the system is performed by 
MD, however, the positional order has oscillations be- 
cause of the momentum conservation. This oscillation 
prevents us from studying the order parameter, there- 
fore, MC simulations are performed in this study. 

Each system contains N particles with the radius a. 
The density is normalized to be p = 1 when the system is 
in the perfect square/cubic lattice, that is, p = {2d/L)'^N 
with the dimension of the system d and the linear size 
of the system L. Throughout this study, the number of 
particles N = 23288 for two- and N = 32000 for three- 
dimensional systems and up to 512 independent samples 
are averaged for each density. The step length is set to 
be CTs = 0.2cr with the radius a J2j|. At the beginning 
of each run, the particles are set up in the perfect or- 
dered configuration, namely, the hexagonal lattice in the 
two-dimensional and the face centered cubic lattice. The 
periodic boundary conditions are taken along the all axes. 
The densities from p — 0.7 to 1.0 are studied. 

The time evolutions of the mean square displacements 
M are shown in Fig. Q] We can see the normal diffu- 
sions starting after some initial relaxations, for example, 
the data of three-dimensional system with p = 0.9 has 
a bend at i ~ 10'^. It implies that the normal diffusion 
started after the positional order are almost destroyed. 
The time evolutions of the positional order parameters 
are shown in Fig. |21 While the positional order is well 
approximated by the second order cumulant in the re- 
gion where the positional order parameter is not so small, 
there are differences especially in the low densities. 

These differences are caused by the higher order cu- 
mulants which are ignored in Eq. Q- The contribution 
from the higher order cumulants can be explained by a 
swapping diffusion process. In particle systems, there 
are two kinds of ways to diffuse; by normal diffusion and 
the swapping. While the normal diffusion destroys the 
positional order parameters as described in Eq. (0, the 
swapping does not. In Fig. |2| the diffusion behavior in 
two-dimensional solid is shown. The density is p = 0.92 
which is high enough than the melting points pm, i.e., 
Pm ~ 0.902 reported by ZoUwcg and Chester |g and 
Pm < 0.905 by Weber et al. |3, and p™ ~ 0.893 by 
Watanabe et al. [l3|- The diffusion shows logarithmic 
behavior up to t ~ lO'' as Mermin predicted Pf- How- 
ever, it varies from the logarithmic behavior because of 
the swapping diffusion around at t '--^ 10^. The distri- 
bution of the displacement u^ at this time is shown in 
Fig. ^ The points around at the center correspond to 
the results of the normal diffusion and the six groups 
around the center group correspond to that of the swap- 
ping diffusion. 

In order to treat the effect of the swapping, we consider 
the system with two types of diffusion, the continuous 
diffusion and the swapping diffusion with a swapping rate 
Er on the lattice with a lattice constant a. The rate £",- 
denotes the probability to jump to the nearest position at 
equilibrium per unit time. The diffusion with swapping 
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FIG. 3: (a) Mean square displacement M of the two- 
dimensional solid. The number of particles A^ = 23288 
and the density p — 0.92. The solid line is Cilni with 
Ci — 1.6 • 10~^. Error bars are smaller than the size of the 
symbols. It shows that the crossover from the normal diffu- 
sion (logarithmic) to the swapping diffusion (power- law), 
(b) The time evolution of the value ln{^ /'i/'), where ^ is the 
positional order parameter with the definition in Eq. Q and 
'5' is the value calculated from diffusion using Eq. JSJ. The 
decimal logarithms are taken for the both axes. The solid and 
dashed lines are drawn for for the guides to the eyes; The solid 
line is C2t^-^^ and the dashed line is Cst with C2 = 4.0-10"^" 
and C3 — 3.5 ■ 10~^. It shows that the exchanging rate in- 
creases as Er ~ t"'^^. 



(u^) is expressed to be, 

(u^)' = (uf) + da^Ert, 



(11) 



with the diffusion without swapping (uf) [23|. In the 
following, the positional order parameter calculated from 
Eq. |SJ| is denoted by '^' in order to distinguish from 
the original definition in Eq. JQ. Using Er, ^' can be 
expressed to be. 



\I'' = exp 



K^u^)" 



2d 
*cxp(-27r^£;ri). 



(12) 



Therefore, Vf' is always smaller than ^P, since Er > 0. 
The contribution of the higher order cumulants is ex- 
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FIG. 4: Distribution of displacements u^ at i = 5 • 10^ of the 
two-dimensional system. A small system with TV = 2900 is 
shown for visibility. The density of the system is p = 0.92. 
The lattice constant is a = 2 and the radius of particles is 
a — 0.89 in this scale. The center (0, 0) corresponds to the 
initial position Ri. The points around at the center corre- 
spond to the normal diffusion and the six small groups around 
the center group correspond to the swapping diffusion. 



pressed to be, 

\n{-9/^') = 2TT^Ert. (13) 

The time evolution of the value In (^/^') is shown in 



Fig. |3fb). It increases as ~ t^"^^, which is faster than 
linear increase. If the exchanging rate Er is constant, the 
value should increase as ~ t. Therefore, the exchanging 
rate Er increases. It implies that the destruction of the 
positional order enhances the swapping of the particles. 

To summarize, we study the dynamics of the positional 
order in the particle systems based on the diffusion pro- 
cesses. We discuss the relation between the positional 
order parameter ^ and the mean square displacement 
M with the cumulant expansion. We find that there are 
two kinds of diffusion processes in particle systems, one 
is the normal diffusion and another is the swapping diffu- 
sion which allows particles to diffuse without destroying 
the positional order. These diffusion processes can be un- 
derstand as cumulants of the displacements; the normal 
diffusion is the second order cumulant, and the swapping 
diffusion contributes to the higher orders. This swapping 
diffusion process will play important roles in systems with 
two or more kinds of particles in high density region. The 
presented arguments are very general, and applicable to 
other systems with general pair potentials. Studying dy- 
namic aspects of Alder transitions based on the cumulant 
expansion should be a further issue. 

We thank M. Mori for helpful suggestions. The compu- 
tation was partially carried out using the facilities of the 
Supercomputer Center, Institute for Solid State Physics, 
University of Tokyo. This work was supported by the 
21st COE program, "Frontiers of Computational Sci- 
ence" , Nagoya University and the Grant-in- Aid for Sci- 
entific Research (C), No. 15607003, of Japan Society for 
the Promotion of Science. 



[1] B. J. Alder and T. E. Wainwright, Phys. Rev. 127, 359 

(1962); W. W. Wood and J. D. Jacobson, J. Chem. Phys. 

27, 1207 (1957). 
[2] N. D. Mermin, Phys. Rev. 176, 250 (1968). 
[3] B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. 41, 

121 (1978); A. P. Young, Phys. Rev. B 19, 1855 (1979). 
[4] J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1181 

(1973); J. M. Kosterlitz, J. Phys. C 7, 1046 (1974). 
[5] S. T. Chui, Phys. Rev. Lett. 48, 933 (1982). 
[6] J. A. ZoUweg, G. V. Chester, Phys. Rev. B 46, 11186 

(1992). 
[7] J. Lee and K. J. Strandburg, Phys. Rev. B 46, 11190 

(1992). 
[8] H. Weber, D. Marx and K. Binder, Phys. Rev. B 51, 

14636 (1995). 
[9] J. F. Fernandez, J. J. Alonso and J. Stankiewicz, Phys. 

Rev. Lett. 75, 3477 (1995). 
[10] A. Jaster, Phys. Rev. E 59, 2594 (1999). 
[11] S. Sengupta, P. Nielaba, and K. Binder, Phys. Rev. E 

61, 6294 (2000). 
[12] K. J. Strandburg, Rev. Mod. Phys. 60, 161 (1988). 
[13] H. Watanabe, S. Yukawa, Y. Ozeki, and N. Ito, Phys. 

Rev. E 66, 041110 (2002); H. Watanabe, S. Yukawa, Y. 

Ozeki, and N. Ito, Phys. Rev. E 69, 045103(R) (2004); 



H. Watanabe, S. Yukawa, and N. Ito, Phys. Rev. E 71, 
016702 (2005). 

[14] N. Ito, Pramana-J. Phys. 64, 871 (2005). 

[15] K. Zahn and G. Maret, Phys. Rev. Lett. 85, 3656 (2000). 

[16] In the low density limit, the mean square displacement 
M of MG is exactly obtained to be M = a^d/{d + 2) ■ t 
with the step length Ua and the dimension of the system 
d. Since the step length is set to be CTs = 0.2a in this 
study, the lines 0.04d/(d + 2) ■ t are shown in Fig.0 

[17] B. Jancovici, Phys. Rev. Lett. 19, 20 (1967). 

[18] D. C. Rapaport, J. Gomput. Phys. 34, 184 (1980). 

[19] B. D. Lubachevsky, J. Gomput. Phys. 94, 255 (1991). 

[20] M. Isobe, Int. J. Mod. Phys. G 10, 1281 (1999). 

[21] We have also performed simulations with other values of 
step length CTs — 0.01,0.05 and 0.1, and confirmed that 
the parameter only changes the time scale of the system 
and the results are not changed qualitatively. 

[22] Consider a diffusion process on the grid with the lattice 
constant a. If particles jump randomly to the next site 
every r, the system involves diffusion with a diffusion 
constant D — a^/2r and the exchanging rate Er is de- 
noted by Er — T^^ . 



